Noname manuscript No. 

(will be inserted by the editor) 



Stability of the splay state in networks of 
pulse-coupled neurons 

Simona Olmi • Antonio Politi • Alessandro 
Torcini 



the date of receipt and acceptance should be inserted later 

Abstract We analytically investigate the stability of splay states in networks 
of iV pulse-coupled phase-like models of neurons. By developing a perturbative 
technique, we find that, in the limit of large N, the Floquet spectrum scales as 
l/N 2 for generic discontinuous velocity fields. Moreover, the stability of the so- 
called short-wavelength component is determined by the sign of the jump at the 
discontinuity. Altogether, the form of the spectrum depends on the pulse shape 
but is independent of the velocity field. 
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1 Introduction 

The first objective of a (neural) network theory is the identification of the asymp- 
totic regimes. The last-decades activity have led to the discovery of fully- and 
partially-synchronized states, clusters and splay or asynchronous states in pulse- 
coupled networks [1,2,3,4 . It has also been made clear that ingredients such as 
disorder (diversity of the neurons and structure of the connections) are very im- 
portant in determining the asymptotic behaviour, as well as the possible presence 
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of delayed interactions and plasticity J5][6]. However, even if one restricts the anal- 
ysis to identical, globally-coupled oscillators, there are very few theoretical results: 
they mostly concern fully synchronized regime or specific types of neurons (e.g. 
the leaky integrate-and-fire model) [U[7|[8]. 

In this paper, we develop a perturbative analysis for the stability of splay 
states (also known as antiphase states [5], "ponies on a merry-go-round" |l(Jj . or 
rotating waves in ensembles of N identical fully pulse-coupled neurons. In a 
splay state, all the neurons follow the same periodic dynamics except for a time 
shift that is evenly distributed. Splay states have been identified in experimental 
measurements performed on electronic circuits [11] and on multimode lasers |12j . 
Theoretical studies have been devoted to splay states in fully coupled Ginzburg- 
Landau equations [13], Josephson arrays [14] . laser models [15], traffic models [16] . 
and pulse-coupled neuronal networks [2]. In the latter context, splay states have 
been mainly investigated in leaky-integrate-and-fire (LIF) neurons [2,3, 17, 18 , but 
some studies have been also devoted to the 0-neurons [19j and to more realistic 
neuronal models [20] . Finally, splay states are important in that they provide the 
simplest instance of asynchronous behaviour and can be thereby used as a testing 
ground for the stability of a more general class of dynamical regimes. 

Our model neurons are characterized by a membrane potential u that is contin- 
uously driven by the velocity field F(u), from the resetting value u = towards the 
threshold u = 1 (see the next section for a more precise definition). As threshold 
and resetting value can be identified with one another and thereby u interpreted 
as a phase, it will be customary to refer to the case F(l) / ^(0) as to that of a 
discontinuous velocity field. Additionally, we assume that the post-synaptic poten- 
tial (PSP) has a stereotyped shape, the so-called a-pulse, that is characterized by 
an identical rise and decay time l/a [2] . As already discussed in [18] , the Floquet 
spectrum is composed of two components: (i) long wavelengths (LWs), which can 
be studied in terms of a suitable functional equation for the probability distribu- 
tion of the membrane potential u [2]; (ii) short-wavelengths (SWs), which typically 
correspond to marginally stable directions in the thermodynamic limit (N — » oo). 
By developing an approach that is valid for arbitrary coupling strength and is 
perturbative in the inverse system-size l/N, we prove that the SW component of 
the Floquet spectrum scales as l/N 2 and is proportional to F(l) — F(0), i.e. it is 
present only if the velocity field is discontinuous. We are also able to determine 
the spectral shape and find it to be universal, i.e. independent of the details of the 
velocity field. 

More precisely, we first build the corresponding event-driven map, by expand- 
ing it in powers of l/N (a posteriori, we have verified that it is necessary to reach 
the fourth order). Afterwards, the expression of the splay state is determined: this 
task corresponds to finding a fixed point of the event-driven map in a suitably 
moving reference frame - analogously to what previously done in specific contexts 
18,21,22 . In practice this task is carried out by first taking the continuum limit 
for the various orders and obtaining suitable differential equations, whose solution 
allows proving that all finite-size corrections for both the period T and the mem- 
brane potential vanish up to the third order. Next, the stability analysis is carried 
out to determine the leading term to the Floquet spectrum. This task involves the 
introduction of a suitable Ansatz to decompose each eigenvector into the linear 
superposition of a slow and a rapidly oscillating component. The following contin- 
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uum limit shows that the two components satisfy an ordinary and a differential 
equation, respectively. 

Altogether, the proof of our main result requires determining all terms up to the 
third order in the 1/N expansion of the splay state solution, while some third order 
terms are not necessary for the tangent space analysis. Going beyond discontinuous 
fields would require extending our analysis to account for higher order terms and 
this might not even be sufficient to characterize analytic velocity fields. In fact, 
previous numerical simulations [21] suggest that the Floquet exponents scale with 
higher powers of 1/N that depend on which derivatives of F(u) are eventually 
discontinuous. Moreover, it is worth recalling that in the case of a strictly sinusoidal 
field, a theorem proved by Watanabe and Strogatz 23; implies that N — 3 Floquet 
exponents (N — 2 for a splay solution) vanish exactly for any value of N. 

In the small coupling limit, one can combine our results with those of Abbott 
and van Vreeswijk 2_ (that are valid only in that regime) for the LW spectral 
component and conclude that the splay state is stable whenever F(0) > F(l) and 
the pulses are sufficiently broad, for excitatory coupling, while it is always unstable 
for inhibitory coupling and any finite pulse-width. This scenario is partially rem- 
iniscent of the stability of synchronous and clustered regimes that is determined 
by the sign of the first derivative dF/du of the velocity-field averaged on the in- 
terval [0, 1] (this latter problem has been investigated in excitatory pulse-coupled 
integrate-and-fire oscillators subject to 5-pulses [Tll24p. 

Section II is devoted to the introduction of the model and to a brief presentation 
of the main results, including an expression for the leading correction to the period 
for the LIF model, to provide evidence that they are typically of 4th order. A 
general perturbative expression for the map is derived in Sec. Ill, while Sec. IV is 
devoted to deriving the splay-state solution up to the third order in 1 /N. The main 
result of the paper is discussed in Sect. V, where the Floquet spectra are finally 
obtained. Sect. VI contains some general remarks and a discussion of the open 
problems. The technical details of some lengthy calculations have been confined 
in the appendices: Appendix A is devoted to the derivation of the splay state 
solution; Appendix B contains the derivation of the leading term (of order four) of 
the period T for the LIF model; Appendix C is concerned with the linear stability 
analysis. 

2 Model and main results 

We consider a network of N identical neurons (rotators) coupled via a mean-field 
term. The dynamics of the i-th neuron writes as 

iii(t) = F( Ui ) + gE(t) = Fi(t) i = l,...,N , (1) 

where in(t) represents the membrane potential, E(t) is the forcing field, and g is 
the coupling constant. When the membrane potential reaches the threshold value 
Ui(t) = 1, a spike is sent to all neurons (see below for the relationship between 
the single spikes and the global forcing field E(t)) and it is reset to Ui(t) = 0. The 
resetting procedure is an approximate way to describe the discharge mechanism 
operating in real neurons. The function F represents a velocity field for the isolated 
neuron and it is assumed to be everywhere positive (thus ensuring that the neurons 
repetitively fire, since they are supra-threshold) , while Ti is the velocity field seen 
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by the neuron i in the presence of a coupling with other neurons. While we consider 
both excitatory (g > 0) and inhibitory networks (g < 0), it is easy to show that 
T remains always positive to ensure the existence of splay states. For the simple 
choice 

F(u) = a - u (2) 

the model reduces to the well known case of LIF neurons. 

The field E is the linear superposition of the pulses emitted in the past when 
the membrane potential of each single neuron has reached the threshold value. By 
following Ref. [2 , we assume that the shape of a pulse emitted at time t = is 
given by E s (t) = ^^e - " , where 1/a is the pulse-width. This is equivalent to 
saying that the total field evolves according to the equation 

2 

E(t) + 2aE(t) + a 2 E(t) = — ( 3 ) 

n\t n <t 

where the sum in the r.h.s. represents the source term due to the spikes emitted 
at times t n < t. 

It is convenient to transform the continuous-time model into a discrete-time 
mapping. We do so by integrating the equations of motion from time t n to time 
(where t n is the time immediately after the n-th pulse has been emitted). 
The resulting map for the field variables reads, 

E n+1 = [E n + r„P„]e- QT " , (4) 
a 2 

1 1 i ) —OLT n . LX 

where t„ = t n +i — t n is the interspike time interval and, for the sake of simplicity, 
we have introduced the new variable P := aE + E. 

In this paper we focus on a specific solution of the networks dynamics, namely 
on splay states, which are asynchronous states, where all neurons fire periodically 
with period T and two successive spike emissions occur at regular intervals t„ = 
T/N. The first result of this paper is that under the assumption that the velocity 
field F(u) is differentiable at least four times, the dependence of the period T onto 
the size iV is of order o(l/N ). In the specific case of LIF neurons, we show in 
Appendix B that the leading correction ST to the infinite size result is indeed of 
order 0{1/N 4 ) and, more precisely, 

_ Kja) - 6 [a(l - e- T ) - l] T 5 

720 ge- T + a(T+ 1 - e~ T ) - 1 TV 4 ' 1 ' 

where K(u) encodes the information on the pulse dynamics (see Eq. (|66p l. We 
did not dare to estimate the quartic contribution for generic velocity fields, not 
only because the algebra would be utterly complicated, but also since our main 
motivation is to determine the leading contributions in the stability analysis, and 
it turns out that it is sufficient to determine the splay state up to the third order. 

The study of the stability requires determining the Floquet spectrum, i.e. the 
complex eigenvalues of a given periodic orbit of period T. With reference to a 
system of size N, the Floquet multipliers can be written as 

Mfc = e «*» e (X *" Hi " k)T " , k = 1, . . . ,3JV , (6) 
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where <f>k represents the Oth order phase (that is responsible for the high frequency 
oscillations of the corresponding eigenvector - see Sec. V) , while \k and Wk are the 
real and imaginary parts of the Floquet exponent, respectively. In the following 
we prove that the leading term of the SW component (i.e. for <j>k away from zero), 
is 

go? F(l)-F(0) ( 6 \ 1 

k 12 J"(1)-F(0) ^1-cos^fc J iV 2 ' K 1 

For discontinuous velocity fields, the real parts of the spectrum scale as 1/N 2 , 
while the imaginary parts are of even higher order. 

For continuous fields, it has been numerically observed that the scaling of the 
spectrum is at least 0(1/N 4 ) [21]. In other words the shape of the spectrum is uni- 
versal, apart from a multiplicative factor that vanishes if and only if F(l) = F(0), 
i.e. for true phase rotators where u = coincides with it = 1. The stability of the 
splay state can be inferred by the sign of F(1)—F(0): in the case of excitatory (resp. 
inhibitory) coupling, the state is stable whenever F(0) > F(l) (resp. F(0) < F(l)). 
In the limit <f>k — > the expression reported in parenthesis in Eq. (J7]) diverges, in- 
dicating that the perturbative analysis breaks down. This limit corresponds to 
the LW component, where our approach can be complemented by that of Abbott 
and van Vreeswijk 2 , which reveals that the corresponding Floquet exponents do 
not depend on the system size. For sufficiently small couplings (\g\ << 1), they 
also found a condition similar to the one reported above, namely that, irrespec- 
tively of the sign of the coupling, the splay state is stable whenever F(0) > F(l) 
for sufficiently broad pulses. In fact, above a critical a-value (i.e. below a given 
pulsewidth), the splay state looses stability due to a supercritical Hopf bifurca- 
tion, which leads to the emergence of a more complex collective regime, termed 
partial synchronization [3ll25| . By combining the conditions for the SW and the 
LW spectrum, one can predict the overall stability of the splay state. In particu- 
lar, the state is stable for excitatory coupling if F(0) > F(l) (and a sufficiently 
small), while it is always unstable for finite networks, for inhibitory coupling, since 
the SW and LW stability conditions are opposite to one another. This last result 
is consistent with the findings reported by van Vreeswijk for inhibitory coupling 
and finite pulse width [3]. 



3 Event driven map 

By following Ref. 26,21 , it is convenient to pass from a continuous to a discrete 
time evolution rule, by introducing the event-driven map which connects the net- 
work configuration at subsequent spike emissions occurring at time t n and t n +i- 
The membrane-potential value Ui{t~ +l ) just before the emission of the (n + l)-th 
spike can be obtained by formally integrating Eq. ([TJ) , 

dtF{ Ul (t))+g / dt[E n +P n {t-t n )]e- at = A1+A2 , 

(8) 

where the minus superscript means that the map construction has not yet been 
completed. This task is accomplished by ordering the membrane potentials from 
the largest (j = 1) to the smallest value (j = N) value and by passing to a 
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comoving frame that advances with the firing neuron, i.e. by shifting the neuron 
index by one unit, 

u n+1 j-i = uj(t~ +1 ) , (9) 

where the first subscript indicates that the variable is determined at time t n +i- 
This change of reference frame allows treating the splay state as a fixed point of 
the event driven map. 

The first integral appearing on the rhs of Eq. (JSj) is now solved perturbatively 
by introducing a polynomial expansion of Ui(t) around t = t n , which, up to third 
order, reads as 



j(t) = u n ,j + iinjSt + ^u n jSt 2 + ^u n ,j5t 3 + O (5t 4 ) 



(10) 



where St = t — t n . Explicit expressions for the time derivates of Uj can be obtained 
from Eq. ([TJ) and its time derivatives, 

U n ,j = F (Un,j) Un, 3 + gE n , 
U n ,j = F" (U n ,j)u n j + F' {U n ,j)u n ,j + gE n , 



where one can further eliminate E n with the help of Eq. ([3} . 

By inserting the expansion (|10[) into the expression of Ai , expanding the func- 
tion F(u), and performing the trivial integrations, one obtains 



r 2 j- 

Al = Fn^jTn ~\~ F n jJ~ n,j~^~ ~\~ ^ 



1 Tn 

Fn,j + gEnF^jJ + {F"[jT, 



■3 



(11) 



,jFn,j*F n> j + F n jJ- n j + g ^^3F n jJ- n j + F n j aF n j^j E n aF nt 

where T n = tn+i — t n and we have introduced the short-hand notation F n ,j for 
F(u n ,j) (and analogously for J 7 ). 

The explicit expression of Ai reads 



a a or 



(12) 



gE n r n + gE n ^ - ga( E n + P n ) -f + ga z ( E n + 2P n ) ■£ + 0(r°) 



21 



Now, by assembling Eqs. (|8l9lllll2p . we obtain the final expression for the evolu- 
tion rule of the membrane potential, 



«n+l,j-l — "nj + FnjTn + [gE n + F^jJ^nJ^ + {F„j [F n jF n J + 



+KjJ%j - 9<x[P n + E n ] } ^ + {-ga(E n + P n )F nij + iK,iKi^ 



~\~F n jJ- nj + g 



(13) 



SF^Fnj + K,i) En + a (E„ + 2P„)j + KjJ%j } ^ + 0{t° 



Eqs. ^ and (|13p define the map we are going to investigate in the following 
sections. The time needed to reach the threshold r n can be determined implicitely 
from Eq. (|13p by setting j = 1, since by definition of the model u nj o = 1. 
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4 Splay state solution 

The splay state is a fixed point of the previous mapping corresponding to a constant 
interspike interval r = T/N. Since the fixed point solutions do not depend on the 
index n they are denoted as, 

E n = E , P n = P , u n ,j = Uj . (14) 

In order to study the dependence of the splay state on the system size N, it 
is necessary and sufficient to formally expand the expression of the membrane 
potentials as follows 

h=0,4 

and, analogously, for the period T, 



yW / 1 \ 

h=OA 



This expansion can be performed by exploiting the explicit dependence of E and 
P on N, as detailed in Eqs. (|56l55p in Appendix lAl 

Finally by substituting the expressions (|16I55I56I57|) in Eq. (|13[) one obtains 
the evolution equations for the membrane potentials 



Ah) _ -(h) {h) 
h=0A h=lA 



N h ^ N h \N 5 J ' 



where the Q variables are defined in Appendix |A"1 

In the large N limit, one can introduce the continuous spatial coordinate x = 
j/N. In practice, this is tantamount to write, 

U {h \x = j/N) = u ( j h) , h = (),•••, 4 . (18) 

It is important to stress that the event-driven neuronal evolution in the comoving 
frame implies that U(0) = 1, i.e. the first neuron will fire at the next step, and 
1/(1) = 0, i.e. the membrane potential of the last neuron has been just reset to 
zero. This implies that C/ (o) (0) = 1 and t/ (0) (l) = 0, while U {h} {0) = U {h) {l) = 
for any h > 0. 

Furthermore, by expanding U^ h '(x) around x = j/N, one obtains 

= -!/*) = *»(.)+ £ -L^)" + o(^), 

m=l,4 

(19) 

By inserting this expansion into Eq. (|17p . we obtain an equation that can be 
effectively split into terms of different order that will be analysed separately. Notice 
that by retaining terms of order h, it is possible to determine the original variables 
at order h — 1. 
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4.1 Zeroth order approximation 

By assembling the first order terms, we obtain the evolution equation for the zeroth 
order membrane potential, namely 

^ = _ 5 _ T (0) F(C /(0) ) . ( 20) 

This equation is equal to the evolution equation of the membrane potential for a 
constant field E, with x playing the role of a (inverse) time. Please notice that, 
up to first order, E = l/T^ ' (see Eq. (|56p ). An implicit and formal solution of 
Eq. (EDJ) is, 

where we have imposed the condition £A°^(1) = 0. However, there is a second 
condition to impose, namely U^°\o) = 1. This second condition transforms itself 
in the equation defining the interspike time interval T^°\ when N — > oo (i.e. in 
the thermodynamic limit) 

dU {0) 

(22) 



o g + T(°)F(UW) 



This result is, so far, quite standard and could have been easily obtained by just 
assuming a constant field E in equation (QJ. If we introduce the formal relation 

F'[U^{x)] = ^ffi in E{ 1- (EHD we obtain 

g+ V*Z°>r - F ' [um(x)]dx • (23) 

which can be easily integrated 



(24) 



giving the following relation (already derived in [25] , by following a different ap- 
proach) 

e -T<°>ff(0) e -T'°»H(l) 

T{U{Q)) = T{U(1)) ' (25) 
where, for later convenience, we have introduced 

H(x)= [ X F'[U {0 \y)]dy , (26) 



and where, for the sake of simplicity, the prime denotes derivative with respect to 
the variable and the dependence of F and F' on U ^ has been dropped. 
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4.2 First order approximation 

By collecting the terms of order 1/N 2 , one obtains 



du(1) = + 1 d2lj{0) ft (1) - i(T {0) ) 2 F'F-^-T {0) F' (27) 

dx 2 dx 2 2 y ' 2 ' v ' 

An explicit expression for the second derivative of (x) appearing in Eq. (|27|) 
can be computed by deriving Eq. (|20[) with respect to x. This allows rewriting 
Eq. (|27[) in a simplified form, namely 

dU^jx) = _ u (i) T (0) F > _ T W F _ (2g) 
dx 

By imposing £/ (1) (l) = 0, one obtains the general solution of Eq. (|28[l , 



U {1) {x)= I duT {1) F[U {0 \u)}exp T W (h{x) - H(u) 

J X 



(29) 



where H(x) is defined by Eq. ()26[) . The further condition to be satisfied, £7^(0) = 
0, implies T^ 1 ' = and thereby we have U^'(x) = 0, i.e. first-order corrections 
vanish both for the period and the membrane potential. 



4.3 Second order approximation 

The second order corrections can be estimated by assembling terms of order 1 /N 3 
and by imposing the previously determined conditions T^ 1 ' = and U^'(x) = 0, 

dU(2) = _ T (0) F ' [/ (2) _ 1 d 3 U {0) _ ^(2) _ (f_ T (Q) F „ 

dx 6 dx 3 6 

_£(T (0) ) 2 [2FF" + F' 2 ] - ^ T( ° ) ' )3 [F"F 2 + F' 2 F] 

Once evaluated d 3 J7 (0) /dx 3 from Eq. ([20]), the above ODE reduces to 

^L = _ u m T m F/ _ T {2) F (30) 

dx 

which has the same structure as Eq. (|28p . Since one has also to impose the same 
boundary conditions as for the first order, namely EA 2 ^(Q) = U^ 2 \l) = 0, we can 
conclude that = and, consequently, U^'ix) = 0. Therefore, second order 
corrections are absent too. 
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4.4 Third order approximation 

By assembling terms of order 1/N 4 , once imposed that first and second order 
corrections vanish, one obtains 

dU(Z) = _ T (0)p'rr(3) , 1 d 4 U m _ (3) _ ^_ (0) ,„ _ ff_ (0) 2 , „ 

dx + 24 di 4 24 & y ' 



- — {T^) 2 FF'" - -L( T <°>) 3 [F' 3 + 8FF'F" + 3gF 2 F"'] 



/ T (0)\4 

24 FF> ^ + AFF " + F3F "'] ■ ( 31 ) 



\4 



By replacing d A U^ /dx 4 with its expression derived from Eq. (|20[) . equa- 
tion (|3ip takes the same form as in the two previous examined cases, namely 

dU{3) _[/(3) T (0) F '_ T (3) J? _ (32) 



dx 

Therefore, we can safely conclude that third order terms vanish too. 

The LIF model can be solved exactly for any value of N, starting from the 
asymptotic value (N — > oo). As shown in Appendix B, it turns out that the 
leading corrections are of fourth order for both the period T and the membrane 
potential. 



5 Linear stability analysis 

The fixed-point analysis has revealed that the finite-size corrections to the station- 
ary solutions are of order o(l/N 3 ). Since such deviations do not affect the leading 
terms of the linear stability analysis (as it can be verified a posteriori) they will 
be simply neglected. Therefore, for the sake of simplicity, from now on, and 
Uj will be simply referred to as T and Uj. 

The evolution rule in tangent space is obtained by differentiating Eq. (|13l) and 
Eq. Q around the fixed point solution. The explicit expression of the correspond- 
ing event-driven map is reported in Appendix C. It consists of evolution equations 
for SP n and 5E n (Eqs. fl77J) and ([78])), and for 5u n ,j (Eq. fl79j)). Finally, 5r n is 
determined from Eq. (|80[) . 

As usual, the eigenvalue problem can be solved by introducing the Ansatz, 

5u n , 3 = fJ-kSuj SP n = HkSP SE n = tik&E 5r n = ^St , (33) 
where /ifc labels the eigenvalues, which must also be expanded as, 

Mfc = e ** e (*-H->)T/* = ^ L + £ £^ + o ( Jj) V (34) 

V h = l,3 ) 

where T^ h > is, in principle, a complex number and, for the sake of simplicity, we 
have dropped its dependence on k. Finally, as already shown, at zeroth order, the 
eigenvalues correspond to a pure rotation (specified by (f>k) with no expansion or 
contraction, i.e. _T*- ^ = 0. 



Stability of the splay state in networks of pulse-coupled neurons 



11 



By inserting the above Ansatze in the map expression (|77I78I79I80|) . one ob- 
tains, after eliminating 5P, 5E and St, a closed equation for Suj, 



1 + 



r {2) r {3) 



N 



+ 



N 2 + N 3 



5u 



+ 



6N 3 J 



5m 



T 



2N 2 



(35) 



T"y 2 f' 2 - 

j 3 J 3 , 3 

2 + 2 ^ 



g_ 2 e 2l4>k + lOe 1 ^ 
+ T° 12(e*** 



1)2 

5 a 2 n(1) 2e 2 ^ -3e^ 



+ 1 



+ 



-T 



+ - — r 

t- 3T 2 



+ 



5e^ + 1 



l) 3 



D--l]+~a 



g 2 e 2 ^" + 10e l4 > k + 1 T 3 



T 12(e^-l) 2 



F'^-F' 



+ 



T 12(e^» - l) 2 
ga 3 e^(e 10 + 1 



T (e^-1) 3 



Su i 



that is the object of our investigation. The overline means that the function is 
evaluated in u^, corresponding to the infinite TV limit. 



5.1 Continuum limit 

Similarly to the splay-state estimation, it is convenient to take the continuum limit. 
However, at variance with the previous case, now one should take in to account 
also the presence of fast scales associated to the "spatial" dependence of 

Therefore, the correct Ansatz is slightly more complicated and we have to 
separate slowly and rapidly oscillating terms, 

S Uj = TTj + ^ 3 e i<l>kJ 1 (36) 

where the complex exponential term accounts for the fast oscillations of the eigen- 
vectors, while, 




h=0,3 h=0,3 



are slowly varying variables. 

Now, we can finally introduce the continuous variable x = j/N, as previously 
done in real space (see Eq. (|18|)). 

n?\x = jj)=*? l) , ef\ x = t) = 'df , (38) 

where h = 0, • • • ,3. This allows expanding bu d -\ = 713—1 + , &j—ie % ^ k ^~ 1 ' around 
x = j/N, similarly to what done in Eq. (|19|) . At variance with the computation 
of the fixed point, now there are also terms like U(l/N) and 6U(1/N), whose 
computation requires a similar expansion but around x = 0. By incorporating all 
the expansion terms within Eq. (|35[) . we have finally an equation, where terms 
of different orders are naturally separated from one another. The calculations are 
summarized in Appendix C and the final equation is (|85[) . By separately treating 
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the different orders, we obtain differential and ordinary equations for the and 77 
variables. It turns out that it is necessary to consider in parallel different orders in 
the fast and slow terms to obtain and 77 to the same order. As a consequence, 
we will see that it is sufficient to expand 5U(1/N) up to order 0(1/N ). 



5.2 Zeroth order approximation 

By assembling terms of order 0(l/N) in Eq. (|85[l . multiplied by the fast oscillating 
factor e 1 ^ 1 -', we obtain a first-order linear differential equation for 0^°\ namely 

d ° (0) ^VfTV'fTTf^- rW 



dx 



m {TF'(U(x))-r w ), (39) 



where T^ 1 -* is the first order correction to the Floquet exponent which should 
be determined. It is important to remind that the prime denotes derivative with 
respect to the variable L^ ', which has been simply redifined U, as previously 
mentioned. The solution is 

{o \x) = K {0) exp \r (1) x - TH(x)] , (40) 

where we made use of the definition (|26p and is a suitable integration con- 

stant. 

By assembling now the slow terms of zeroth order and reminding the definition 
of T(U(x)), we find the following algebric equation 

77W(x)(e^-l) = -[ e ^°)(0) + 77^(0)] y§M ■ (41) 



With the help of Eq. (ji0|) , we obtain 

n W {Q) = _ O i0) {Q) = _ K W e -TH { 0) ^ 
rr(0)/^ _ jAO) -TH{0) F{U{x)) 

( x >~ K e F(U(0)) 

We can now impose the boundary condition 5U {0) (x = 1) = 6> (0) (1) + 77 (0) (1) = 0. 
This implies that 

e -T//(i)+r (1 > e -TH(0) 

.F(t/(1)) = ^([/(0)) ' (42) 

By now exploiting Eq. ()25[) . we find that r^ 1 ' = 0, i.e. the Floquet exponent (both 
its real and its imaginary part) is equal to zero at first order in 1/N. Furthermore, 
Eq. (g0]) becomes 

0^(x) = K^e- TH ^, (43) 

i.e. the eigenvectors are independent of the phase <j>k and are thereby equal to 
one another. In other words we are confirmed that the degeneracy has not been 
removed. 



Stability of the splay state in networks of pulse-coupled neurons 



13 



5.3 First order approximation 

By assembling the fast terms of order 1/iV 2 and by setting = 0, we find that 
0(1) satisfies the following first order differential equation, 



dx 

whose solution is 



' l( ' yl ' r«eM -0^TF'(U(x)) , (44) 



V(x) = (r^K^x + A«) e~ TH ^ , (45) 



where K^ 1 ' is an integration constant associated with the solution of the previous 
equation. 

By collecting the slow terms of order 1 /N in Eq. (|85[l , one obtains the algebric 
equation 



e^0 (1) (O) + 77 (1) (O) 



HU{x)) 



Hum 



(46) 



n {1) (x)(e^ -1) = 
whose solution is, 

77( 1 )(0) = -e (1) (0) = -A"( 1) e- Tff W , 

1 j_ 6 Hum 

By imposing the boundary condition 5U {1 \x = 1) = 6> (1) (1) + 77 (1) (1) = 0, it is 
possible to evaluate r^ 2 \ 

e W (1 ) + 77( 1 )(1) = (JC^r^ + K^)e- TH ^ - K^e- TH ^W^ = . 

(47) 

By again exploiting Eq. (f25|) . we find that r^ 2 ' = and, thereby (from Eq. (US}) 

0( 1) (x) = A-( 1 ) e - TH ^ . (48) 

Altogether, we can conclude that the second order correction to the Floquet ex- 
ponent vanishes as well, and one cannot remove the degeneracy among the eigen- 
vectors. 



5.4 Second order approximation 

By assembling fast terms of order 1/JV 3 appearing in Eq. t)85|) and by setting 
= = 0, the following first order differential equation for can be 

derived 



d ° {2) "(3)0(0) , «(2W 



dx 

whose solution is 



r w G m + 0^'TF'(U(x)) , (49) 



9 {2 \x) = (r^K^x + K {2 ^ e~ THM , (50) 
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where K^ 2 ' is an integration constant associated with the solution of the previous 
differential equation. 

Furthermore, by collecting the slow terms of order 1/N 2 , we obtain the algebric 
equation, 

r (»)w^ e ^ + iOe^ + l F(U(0))-F(U(x)) 



n yz >(x)(e lipk - 1) = ga z T0 {v >(O) 



12(e 1 ^ - 1) [T(U(0))]< 



e i0fc <9 (2) (O) + i7 (2) (O) 



Hum ' 

By imposing that the above equation is satisfied for x = 0, it reduces to 

7j( 2 )(o) = -e( 2 )(o) = -^ (2) e - TH(0) , 

(x)-galU [V) 12(e i^_ 1)2 [^(f/(0))] 2 U W .F(t/(0)) 

Finally, by imposing the boundary condition 5U {2) {x = 1) = 6> (2) (1)+77 (2) (1) = 
0, it is possible to determine J^ 3 - 1 , 

r (3) = ga 2 T F(U(0))-F(U(l)) ( 6 _ A 



12 _F([/(0)).F(E/(1)) \l-cos<j) k 

Accordingly, is real and depends on the difference between F(U(x = 1)) = 
F(0) and F(U(x = 0)) = ^(1), confirming the numerical findings in [21] . There- 
fore, the imaginary terms uii are smaller than 1/N 2 . 

In the specific example of a leaky integrate-and-fire neuron the expression for 
reduces to 



r 



(3) 



^ r ( J -« r -- T )(r^fc-0 ■ < 52 > 



since, by using the equations that characterize LIF neurons, the following relation 
holds 

F(U{1))-F(U(0)) _ 1 _ 1 = (e T +e -T_ 2) m 

All in all, Eq. (|5ip generalizes the expression found for the LIF model Eq. ([52]) _2lJ3- 



6 Conclusions 

We have derived analytically the short- wavelength component of the Floquet spec- 
trum of the splay solution in a finite, fully coupled, network composed of generic 
suprathreshold pulse-coupled phase-like neurons. This component is marginally 
stable in the thermodynamic limit and thereby requires a particular care. The 
analytical estimation of the long-wavelength component was previously derived in 
the small-coupling limit [2] . It would be nice to extend such analysis to finite cou- 
pling strength, but this is a rather problematic goal, since the eigenvalues remain 

1 In comparing with 1211 one should pay attention to the different normalization used here 
to define fi^ in Eq. J34D 
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finite in the thermodynamic limit and so there is no evident smallness parameter 
to invoke for a safe expansion. 

Our analysis has revealed that, in discontinuous velocity fields, the SW spec- 
trum scales as l/N 2 , and the stability is controlled by the sign of the difference 
between the velocity at reset and at threshold. The shape of the spectrum is oth- 
erwise universal, at least for a given choice of the post-synaptic potential. Our 
formalism could be easily implemented for any pulse shape, provided that Eq. ([3]) 
is replaced by the appropriate evolution equation. Preliminary numerical stud- 
ies anyway suggest that different (e.g., purely exponential) pulses yield the same 
scaling behaviour, but are characterized by different Floquet spectra [27] . 

Moreover it is worth recalling that <5-like pulses in networks of LIF neurons give 
rise to a different scenario, with a finite (in)stability of the whole SW component 
[18] . The difference is so strong that the two scenarios cannot be reconciled even 
by taking the limit a — > oo (zero pulsewidth) as the limits N — > oo and zero pulse- 
width limit do not commute [TH]. This reveals that even the simple construction 
of a general stability theory of the splay states requires some further progress. 
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A Fixed-point expansion (general case) 

A simple calculation shows that the splay state expression (114 I t can be easily obtained by first 
solving Eq. [4] (see also 1181 ) 



1 



N (i. 



-aT/N ) 



N (e<*T/N _ i) 



(54) 



where T is the period of the splay state, which must be determined self-consistently. 
The 1/JV expansion of these exact expressions leads to 



a 
T<°) 



'a rW " 


a 


2 T(°) 2 


N + _ 



Q 2 T (o) T (2) T ( X y 



12 



T(°) T(°) z 



a 
AT2 



X 2 T (1) T (3) 



T (i) T (2) T (iy 



12 



T(.°y 



T(°) 



• "(-) 



N 3 



(55) 



1 

T(°) 



T(°) 2 Af 



j2 T (0) T (2) T (\y 



12 



T(»)' 



T (oy 



1 



o? m T( 3 ) tWt( 2 > 
- — T {1 > — - + 2- 



12 



T (iy 
T(°) z 



1 /I 



(56) 



q2_ a 2 T(°) 
2N + 



6A^ 2 



a 3 TW 
W 3 



(57) 
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where we have reported also the expansion of E that is necessary to pass from expression (1 1 3 1> 
to H17I I. Please notice that while the membrane potentials and the period are expanded up to 
0(1/N 4 ), as in II15H and J 161) . here we limit the expansion to 0(1/7V 3 ) terms, since the field 
variables appearing in the event-driven map are integrated over an interspike-interval (see pj)l . 

To proceed further, we need also to introduce the expansions of the velocity field and of 
its derivatives, 



j 

N 
i N 



(2) 
j 



a 



(3) 
3 



3 N 2 
~(2) 

3 N 2 



F'(u ] ) = F 



3 TV 3 
F 



2N 2 



3 2N 2 



-(1) - 2 
' 3 3 

iV 3 



JV3 



6^3 



Viv 4 



where the overline means that the function is computed in «j , which corresponds to the 
infinite N limit. 

By replacing the membrane potentials, the period, the self-consistent fields and the velocity 
field with their expansions, the event-driven map (1 1 3 P can be formally rewritten for the splay 
state as (1171) with the introduction of the following auxiliary variables 



Q 



CO 



- T^Fj 



(58) 



Q(2) 



: T^Fj 



9 £j t (0) 



F'. T (0) 



(59) 



Q (3) 



3 3 2 3 2 3 



F 3 F 3 



6 3 



^3 



\ t(°) 2 f/ 


T (0) + 


/ 6 





T^'Fj +T (2) Fj 



2 
(60) 



1(4) 



T (0) 



(sf^ + f.f;" 



^3 



R 3 3 a 3 3 



T (0)-" 
24 



(61) 

T (0) 2 

T (0) 4 
24 



?. fi (.') + ^[ S W]j +J £L fi w + ^' 

J j 2 J 2 J 6 J 

Ft / — //— — < 



T (o) T (i) 



T (1)2 + 2T (0) T (2) 



Fj + FjTW . 



B Fixed-point expansion (LIF model) 

In the case of the LIF neuron (see Eq. Q), the fixed point of the event-driven map reads 

= e~ T Ui + x , (62) 
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where 



Its solution is 



e _T — e~ aT / 

a(l-e- T )+g (e 

a — I \ 



a - I 

NT+jr 



-9 



a - 1 



x- 



(63) 



(64) 



By expanding Eq. II64II for j = and for a generic j, one can derive perturbative expressions 
for the period T and the membrane potential, respectively. Let us start by substituting the 
expressions H161 1551561 1576 in Eqs. (I63D . This leads to the expansion 



M) 



r 4 
24 



r g t 
120 + T 720 



K(a) + 0(1/N 4 ) 



where 



K(a) 



360a ft 



722a 5 + 363a 4 + 5a 2 - 12a 



(65) 



(66) 



(a - 1)2 

accounts for the dependence on the field dynamics. Now, with the help of Eqs. (116165ft and 
expanding the exponential terms up to the fourth order, we obtain a closed equation for the 
interspike interval, 



m = 1 
1 

+ JV3 



1 
iV 1 



T (4)£( T (o)) + ^( T (0)) + T^w[ A) 



hT (2) W (4) +T (3 )>v (4) 



(67) 



where 



C (T(°)) = -(l-e- T 

^ (0) )=(«+4r 



9 T (0) 3 (1 _ K 



120 

_ T (0) 



(l-e- T(0) ) 



while Wp' identifies a term of order l/i\P that is multiplied by TW. Since, while proceeding 
from lower to higher-order terms, we find that TM = (for i < 4), it is not necessary to give 
the explicit expression of the W> functions as they do not contribute at all. 
One can equivalently expand Uj 



(1) (2) (3) (4) 

(0) U j U j U j U 3 ( 9 \ ( 



T(°)ri--ii 



1 

1 



T (2) ?(T (0) )+T (1) 2 (2) 



1 

JV3 



f(T (0) )T (4) +a (T(0) ) + T m z <A) +T (2) Z W +T (3) Z W 



— f(T(°)) 
TV v ; 



(68) 



where 



, (T (°)) = ^(1 - e ^ (0) (^-i))T(°) 3 (f - l) 



p T(°)(i-l) 



(7) 

while we do not provide explicit expressions for Z\ as they turn out to be irrelevant 
Now we are in the position to analyse the different orders. 
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B.l Zeroth Order 

By assembling the terms of order 1 in Eq. (1676 . we obtain 

(. + 5 fe)(.-.— >-.. 

This is an implicit defintion of the asymptotic interspike time T' ' 

T (0) = In — — 

\T(»)(«-l)+s 



(69) 



(70) 



Analogously, we can find an explicit equation for the membrane potential by assembling the 
terms of order 1 in Eq. 1681 



-(•+*; 



l-e 



In the thermodinamic limit the solution for m^°' becomes 



T (0) 



1/(0)/ 



which coincides with Eq. (1216 with F = a — C/( ' . 



(71) 



(72) 



B.2 From first to third order 

By separately assembling the terms of order 1/N Z (for i = 1, 2, 3) in Eq. 16711 . we obtain 

TW((TW)=0 , (73) 

which implies that TW = since £ 7^ 0. Moreover, by assembling the terms of order 1/JV 1 in 
Eq. 16811 . we obtain 

U W = f ( T (0) )T W _ ( 74 ) 

which thereby implies that first, second and third order corrections vanish also for the mem- 
brane potential. 



B.3 Fourth Order 

The order which reveals a different scenario is the fourth one. By assembling the terms of order 
1/N 4 in Eq. (J67J we obtain 

T (4) = _C(TW) 

f(Tt°)) ' K ' 

whose explicit expression is reported in Eq. (£5J . By analogously assembling the terms of order 
1/7V 4 in Eq. EHJ, we obtain 



„(4) =c(T (0) )T (4) +(j(T (0) ) 
which becomes, in the thermodynamic limit, 

= - («+ ^) T« e ^ (0) C-D (:c -l)- s 
-fl- e T(0) ( a; - 1 ))T( ) 3 f^^-l"l 

> v 6 y 



(76) 



T (4) 

T(°) 2 



T^tTO^-i) 



120 



Stability of the splay state in networks of pulse-coupled neurons 



19 



C Expansion in tangent space around the fixed point 

C.l Introduction 

The first equations of the tangent map can be determined by differentiating Eq. ifJJ and thereby 
expanding in powers of r (this is equivalent to expanding in powers of X/N, as the dependence 
of t on N would only generate higher order terms), 

, a 2 a 3 o a 4 . \ . ~ / , a 3 , » 4 , \ . 

| 1 -ar+ — r 2 - — t 3 + — r 4 J SP n +P f -a + a 2 r - — r 2 + — r 3 J <5r„ 

(77) 



a o a 



1 - ar H t 15 

2 6 

v3 



2d 



a 2 2 a 3 , a 4 

-Q_E 1 — (XT H T 7" H 7 

2 6 24 



^ 3 

2 


-^ 4 
6 , 


| SP n 


2 2 
T — 


-a 3 r 3 + 
3 


5 4 

Q T 

24 





(78) 



5t„ 



where the dependence of r on n has been dropped, since we are considering a linearization 
around the splay state. 

By further differentiating Eq. H13II around the fixed point solution, one obtains 

<S«m+l,i-l = <5«n,i + (-Fi<5« n ,i + 5<5E n ) T + [-^i FiSun,i + P, (-Fi<5«n,i + 3<5i? n ) + g<5_B n ] — 
+ 1 2gF"T t 8E n + F'"T\ + F~" Uf^ + gF^ <5« n , t + F 2 + gSE n ) + gF"^E n 



-ga 



(<5P„ + SE n )} T — + i ^ 











F'iTi + gE 


T + 



Fi ?\ + (F j) 2 T i + gFtE - ga(E + P) 



-r—/t-=2, — ///— 3 



(Fi) J 7 i + ^F l F i T\ +F i r i + ga(2a - F t )P 



Finally, <5r n can be determined by differentiating Eq. (1 1 3 I t for i = 1 



2 
(79) 



<5£ n r 

<5t„ = 4 gr - g 



1 % 1 h 

-{F 1 +a) + =-gE 

z -r\ 



_ ; ( SPn f _ 9 T 3 

2 T 2 T 



S (* + *) + J,(W + rf) , - J 1-|^ 1 i 

4* jaf Ql^ j=! + gp) + ^jg 2 ^ 3 + i*? 2 ^ (« + F[) - | fl Fi'i 
F^ F \ 

In order to find the Floquet eigenvalue fj,^, one should substitute the Ansatze (133 I t into 
Eqs. H77I78H . This allows to find explicit expressions for SP and SE as a function of /ij., r and 
St, namely 



(80) 



SP : 



ar /j fc + 1 2 2 a 3 T 3 Mfc(/Jfc + 1) 
1 z 7 + a t M k - — — 



2^-1 



2 tw, - l) 3 



St 



SE : 



Mfc - 1 



or (pk + 1) 2 2 2 M 3° 3t3 Mfc(Mfc + 1) 
2 - 1 * 2 ( Mfe - l) 3 



5t 
T 7 



(81) 



(82) 
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where we have introduced the shorthand notation 



10/Ltfc + 1 



12(/x fe - l) 2 



(83) 



By substituting SP and SE, as given by (18 1 I t and II82I I. into Eq. (1806 . we can express St 
directly in terms of 6u\ 



9a 2 2 9a 2 
— M^t 



Fi (/*fc + 5) , (»k + 1) 
1 a- 



12 0l fc - l) 2 {H k - 1)3 



<5ui 



(84) 



where we exploited the equality Ti = -Fi + ^ which follows from the fact that in the thermo- 
dynamic limit E = Tp . 

By inserting the expressions in Eqs. (18 1 1 1821 18411 into Eq. (75) , we find a single equation 
for the eigenvalues and eigenvectors. 



F i + F i F i r + 
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F — 1 F — 
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p ^{Fi-F^Mk 
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F\' 



By now substituting the /i^ expansion (1 34 II and retaining the leading terms, we obtain Eq. 1351 . 



G.2 N -> 00 limit 

Once the continuous variables 1 138 II have been introduced, it is necessary to estimate U(l/N) 
and 8U(l/N), by expanding such variables around zero. By inserting the resulting expansion 
for U(l/N) into the expressions for F± and T\, we obtain, respectively 
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An analogous procedure for SU(1/N) leads to 
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where CW, C' 1 ), C< 2 ) are defined according to the following equatic 
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C (2) = e ^ fc 
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We now expand SU(1/N) up to the order O (■^■J, thus neglecting higher orders, because they 
contribute to the definition of 77 variable and we need terms at least of order O (j^j) (one 
order lower than needed to define 6>). By inserting the Ansatz (136 6 and the previous expansions 
in Eq. (135 D ■ we finally obtain a closed equation for the eigenvalues and eigenvectors, 
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where we have introduced the shorthand notation B in order to characterize a term of order 
O (j^t)) whose explicit expression is not necessary, since it turns out to contribute to the 
definition of the 77 variable, and it is therefore one order beyond what we need. Moreover, 
notice that the terms appearing within round brackets in the rhs of the above equation can 
be shown to be zero, due to exact algebric cancellations that emerge from the solution of the 
equation order by order. Finally, 

AW(U(x)) = TF'(U(x)), A^{U{x)) = ^ {f"(U(x))T(U(x)) + [F'(U(x))] 2 } , 

A<-»\u(x)) = ^ {f"'(U(x))[T(U(x))] 2 + AF\U(x))F"(U(x))F(U(x)) + [F'(U(x))] 3 } , 



(85) 



B (o Hu{x)) = ^m, 



B 



«([7(x)) = [TF'(U(x))+TF'(U(0))] 



HU{x)) 

Hum 



B^ 2 \U{x)) =T 2 



g 2 e 2 ^ + We^k + 1 F(U{0)) - F(U(x)) F"{U(x)) [T{U{x))} 2 



— a 
T 



12(e**» - l) 2 



[T(U(0))] 2 



F(U(0)) 
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